#--------------------------------------------------------------------------------------------------------------------#
#------------------------------------------------------ Lea Portmann ------------------------------------------------#
#------------------------------------------------------ October 2021 ------------------------------------------------#
#--------- What Makes a Successful Candidate? Political Experience and Low-Information Cues in Elections ------------#
#---------------------------------------------------- Conjoint Experiment -------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------#

# ------------------------------------------------------------ Load packages

library(Hmisc)
library(plyr)
library(dplyr)
library(lme4)
library(tidyverse)
library(vtable)
library(xtable)
library(estimatr)
library(tidyr)
library(sandwich)
library(lmtest)
library(texreg)
library(dotwhisker)
library(margins)
library(effects)
library(gridExtra)

# ------------------------------------------------------------ Read data

rm(list=ls())
setwd(".../...")

load(".../.../d_long.RData")

# -------------------------------------------------------------------------- Balance checks

d_long %>% dplyr::select(gender, exp_arm, age, educ) %>% 
  vtable::sumtable(vars = c("gender", "age", "educ"), group="exp_arm", group.test=TRUE, out="latex")

# ------------------------------------------------------------ Models

d_long$exp_arm_f <- as.factor(d_long$exp_arm)

### Linear models

m1 <- lm_robust(dv_choice ~ name_f + profession_f + age_f + policy_f + children_d, data = d_long, clusters = ResponseId, se_type = "stata")
summary(m1)

m2 <- lm_robust(dv_choice ~ profession_f + age_f + policy_f + children_d + name_f*IDEOL_f, data = d_long, clusters = ResponseId, se_type = "stata")
sum2 <- summary(m2)

m3 <- lm_robust(dv_choice ~ exp_arm_f*name_f + exp_arm_f*profession_f + exp_arm_f*age_f + exp_arm_f*policy_f + exp_arm_f*children_d, data = d_long, clusters = ResponseId, se_type = "stata")
summary(m3)

d_long_exp1 <- d_long %>% filter(exp_arm_f=="1")

m4 <- lm_robust(dv_choice ~ name_f + profession_f + age_f + policy_f + children_d + polexp_fc, data = d_long_exp1, clusters = ResponseId, se_type = "stata")
summary(m4)

### Logit models

m3_logit <- glm(dv_choice ~ exp_arm_f*name_f + exp_arm_f*profession_f + exp_arm_f*age_f + exp_arm_f*policy_f + exp_arm_f*children_d, family = binomial, data = d_long)
coef_m3 <- coeftest(m3_logit, cluster = d_long$ResponseId)
summary(m3_logit)

# ------------------------------------------------------------ Tables

### Main model output

texreg(m3, booktabs=TRUE,
       omit.coef = "bfscantno",
       custom.coef.names = c("Intercept", "Competence information: provided", 
                             "Name: non-Swiss", 
                             "Profession: commercial employee", "Profession: primary school teacher", 
                             "Profession: lawyer", "Profession: social worker", "Profession: salesman", 
                             "Age: 48 years", "Age: 61 years", "Policy: security policy", 
                             "Policy: environmental policy", "Policy: economic policy", 
                             "Children: 3", 
                             "Competence information (provided) x Non-Swiss name","Competence information (provided) x Commercial employee", 
                             "Competence information (provided) x Primary school teacher", "Competence information (provided) x Lawyer", 
                             "Competence information (provided) x Social worker", "Competence information (provided) x Salesman", 
                             "Competence information (provided) x Age (48 years)", "Competence information (provided) x 61 years", 
                             "Competence information (provided) x Security policy", "Competence information (provided) x Environmental policy",
                             "Competence information (provided) x Economic policy", 
                             "Competence information (provided) x 3 children"), 
       caption="Effect of candidate characteristics on probability of preferring candidate",
       include.ci = FALSE,
       single.row = TRUE,
       caption.above = TRUE, stars = c(0.001, 0.01, 0.05, 0.1))

### Logit

texreg(m3_logit, booktabs=TRUE,
       omit.coef = "bfscantno",
       custom.coef.names = c("Intercept", "Competence information: provided", 
                             "Name: non-Swiss", 
                             "Profession: commercial employee", "Profession: primary school teacher", 
                             "Profession: lawyer", "Profession: social worker", "Profession: salesman", 
                             "Age: 48 years", "Age: 61 years", "Policy: security policy", 
                             "Policy: environmental policy", "Policy: economic policy", 
                             "Children: 3", 
                             "Competence information (provided) x Non-Swiss name","Competence information (provided) x Commercial employee", 
       "Competence information (provided) x Primary school teacher", "Competence information (provided) x Lawyer", 
       "Competence information (provided) x Social worker", "Competence information (provided) x Salesman", 
                             "Competence information (provided) x Age (48 years)", "Competence information (provided) x 61 years", 
                             "Competence information (provided) x Security policy", "Competence information (provided) x Environmental policy",
                             "Competence information (provided) x Economic policy", 
                             "Competence information (provided) x 3 children"), 
override.se = list(coef_m3[,2]),
override.pvalues = list(coef_m3[,4]), 
       caption="Effect of candidate characteristics on probability of preferring candidate",
single.row = TRUE,
       caption.above = TRUE, stars = c(0.001, 0.01, 0.05, 0.1))

### Type of competence cues

library(texreg)
texreg(m4, booktabs=TRUE,
       omit.coef = "bfscantno",
       custom.coef.names = c("Intercept", 
                             "Name: Non-Swiss", 
                             "Profession: commercial employee", "Profession: primary school teacher", 
                             "Profession: lawyer", "Profession: social worker", "Profession: salesman", 
                             "Age: 48 years", "Age: 61 years", "Policy: security policy", 
                             "Policy: environmental policy", "Policy: economic policy", 
                             "Children: 3", "Competence cue: member of cantonal parliament",
                             "Competence cue: informed", "Competence cue: intelligent", "Competence cue: leadership qualities", 
                              "Competence cue: president of party in municipality of residence"),
       reorder.coef = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 15, 16, 17, 18, 14),
       caption="Effect on probability of preferring candidate, different types of competence cues",
       include.ci = FALSE,
       single.row = TRUE,
       caption.above = TRUE, stars = c(0.001, 0.01, 0.05, 0.1))

### Ideology

texreg(m2, booktabs=TRUE,
       omit.coef = "bfscantno",
       custom.coef.names = c("Intercept", 
                             "Profession: commercial employee", "Profession: primary school teacher", 
                             "Profession: lawyer", "Profession: social worker", "Profession: salesman", 
                             "Age: 48 years", "Age: 61 years", "Policy: security policy", 
                             "Policy: environmental policy", "Policy: economic policy", 
                             "Children: 3", "Name: Non-Swiss", "Ideology: center", "Ideology: left", 
                             "Non-Swiss name x center", "Non-Swiss name x left"),
       include.ci = FALSE,
       single.row = TRUE,
       #reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 1),
       caption="Effect of candidate name on probability of preferring candidate, by ideological position of respondents",
       caption.above = TRUE, stars = c(0.001, 0.01, 0.05, 0.1))

# ------------------------------------------------------------ Coefficient plot

modelFrame <- data.frame(Variable = rownames(summary(m1)$coef),
                          Coefficient = summary(m1)$coef[, 1],
                          SE = summary(m1)$coef[, 2])

# Add reference categories

Variable <- c("age_35", "policy_jugendpolitik", "profession_hausarzt", 
              "children_2", "name_moser")
Coefficient <- c("", "", "", "", "")
SE <- c("", "", "", "", "")

ref_cat.data <- data.frame(Variable, Coefficient, SE)

allModelFrame <- merge(modelFrame, ref_cat.data, all=TRUE)

allModelFrame <- allModelFrame[-c(1), ] #Remove intercept

allModelFrame <- allModelFrame %>% map_df(rev)

allModelFrame$Coefficient <- as.numeric(allModelFrame$Coefficient)
allModelFrame$SE <- as.numeric(allModelFrame$SE)

#Define labels

my_labels = c("age_35" = "Age (baseline: 35 years old)", 
              "age_f48" = "48 years old",  
              "age_f61" = "61 years old", 
              "children_2" = "Children (baseline: 2 children)", 
              "children_d3" = "3 children", 
              "exp_arm_0" = "Competence information (baseline: No competence information)", 
              "exp_arm_f1" = "Competence information", 
              "exp_arm_f1:name_fMarko Kovac" = "Competence information x Name Kovac",
              "name_fNon-Swiss name" =  "Non-Swiss name", 
              "name_moser" = "Name (baseline: Swiss name)", 
              "policy_fSicherheitspolitik" =  "Security policy", 
              "policy_fUmweltpolitik" =  "Environmental policy", 
              "policy_fWirtschaftspolitik" =  "Economic policy", 
              "policy_jugendpolitik" =  "Policy (baseline: Youth policy)", 
              "profession_fKaufmännischer Angestellter" =  "Commercial clerk", 
              "profession_fPrimarlehrer" =  "Primary school teacher", 
              "profession_fRechtsanwalt" =  "Lawyer", 
              "profession_fSozialarbeiter" =  "Social worker",  
              "profession_fVerkäufer" =  "Sales person", 
              "profession_hausarzt" =  "Profession (baseline: Family doctor)")


# Plot function

PlotCoefficients <- function(z) {
  z <- ggplot(allModelFrame)
  z <- z + geom_hline(yintercept = 0, colour = gray(1/2), lty = 5)
  z <- z + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                              ymax = Coefficient + SE*interval1),
                          lwd = 0.9, position = position_dodge(width = 1))
  z <- z + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval11,
                              ymax = Coefficient + SE*interval11),
                          lwd = 1.4, position = position_dodge(width = 1))
  z <- z + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval1,
                               ymax = Coefficient + SE*interval1),
                           lwd = 0.5, position = position_dodge(width = 1)
  )
  z <- z  + theme_bw() + coord_flip()
  z <- z + ggtitle("")
  return (z)}

#Specify 95% confidence intervals

interval1 <- -qnorm((1-0.95)/2)  # 95% multiplier
interval11 <- -qnorm((1-0.90)/2)  # 90% multiplier

#Reorder variables

order_var <- c("children_d3", "children_2", 
               "policy_fWirtschaftspolitik", "policy_fUmweltpolitik", "policy_fSicherheitspolitik", "policy_jugendpolitik", 
               "age_f61", "age_f48", "age_35", 
               "profession_fVerkäufer", "profession_fKaufmännischer Angestellter", 
               "profession_fSozialarbeiter", "profession_fPrimarlehrer", 
               "profession_fRechtsanwalt", "profession_hausarzt", 
               "name_fNon-Swiss name", "name_moser")

allModelFrame$Variable <- factor(allModelFrame$Variable, levels=order_var)


#Define theme

my_theme <- theme(axis.text.y = element_text(size = 12, colour = "black")) + theme(axis.text.x= element_text(size = 12, colour = "black")) +  
  theme(axis.title = element_text(size = 12)) + theme(legend.position="bottom") + theme(legend.title = element_text(size = 12)) + 
  theme(legend.text = element_text(size = 12)) + theme(plot.margin = unit(c(0,1,0,0), "cm")) 

#Coefficient plot
coef_conjoint <- PlotCoefficients(z)
coef_conjoint <- PlotCoefficients(z) + scale_y_continuous("Effect on probability of preferring candidate", limits = c(-0.25, 0.25)) + scale_x_discrete("", labels = my_labels) +
  my_theme + theme(axis.text.y = element_text(face = c("plain", "bold", "plain", "plain", "plain", "bold", "plain", "plain", "bold", "plain", "plain", "plain", "plain", "plain", "bold", "plain", "bold"))) +
  scale_colour_manual("Model", values = "black") + theme(legend.title = element_blank())
coef_conjoint

jpeg("coef_conjoint.jpeg", 
width = 8, height = 6, units = 'in', res = 800)
plot(coef_conjoint)
dev.off()

# ------------------------------------------------------------ Marginal effects and difference

# ---------------- Competence cues

### Marginal effects of name by experimental arm

m3_marg <- summary(margins(m3, at = list(exp_arm_f = c("0","1")), variables = "name_f", data=d_long))

colors_1 <- c("black", "black")

m3_marg$exp_arm_f <- factor(m3_marg$exp_arm_f, 
                  levels = c(0,1), 
                  labels = c("No competence cues", "Competence cues"))

m3_eff <- ggplot(m3_marg, aes(x = exp_arm_f, ymin=lower, ymax = upper,
                                color=factor)) + geom_errorbar(color=c("black", "black"), width=0.3, size=1, 
                                                                   position=position_dodge(.5)) + 
  geom_point(aes(x=exp_arm_f, y=AME, color=factor), color=colors_1, size=3) +
  theme(legend.position="none") +
  coord_flip() +
  labs(x="", y="Effect of non-Swiss name on probability of being preferred", title="A", color = "Candidate name") + 
  theme_bw() +
  scale_x_discrete(limits = rev(levels(m3_marg$exp_arm_f))) + 
  geom_hline(yintercept=0, linetype="dashed", size=1)  +
  #theme(axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
  scale_fill_manual(values=colors_1) + theme(legend.position="right") + 
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14), plot.title = element_text(size=16)); m3_eff

### Difference Z1 -Z0

x <- summary(m3)$coef

modelFrame <- data.frame(Variable = rownames(summary(m3)$coef),
                         Coefficient = summary(m3)$coef[, 1],
                         SE = summary(m3)$coef[, 2], 
                         lower = summary(m3)$coef[, 5],
                         upper = summary(m3)$coef[, 6])

modelFrame <- modelFrame[15, ]


m3_diff <- ggplot(modelFrame, aes(x = Variable, ymin=lower, ymax = upper)) + 
  geom_errorbar(color=c("black"), width=0.2, size=1, 
  position=position_dodge(.5)) + 
  geom_point(aes(x=Variable, y=Coefficient), color=c("black"), size=3) +
  #scale_y_continuous(limits=c(0, 1)) +
  theme(legend.position="none") +
  coord_flip() +
  labs(x="", y="Difference of effects", 
  title="B", color = "Candidate name") + 
  theme_bw() +
  #scale_x_discrete(limits = rev(levels(m1_marg$exp_arm_f))) + 
  geom_hline(yintercept=0, linetype="dashed", size=1)  +
  theme(legend.position="right") + 
  scale_x_discrete("", labels = "                  Non-Swiss name") + 
  theme(axis.text=element_text(size=14),
   axis.title=element_text(size=14), plot.title = element_text(size=16)); m3_diff

### Figure

m3_eff_diff <- grid.arrange(m3_eff, m3_diff, 
                             ncol=2, nrow = 1); m3_eff_diff 

# ---------------- Ideological position 

### Predicted values

ef <- effect("name_f*IDEOL_f", m2)
summary(ef)
ef <- as.data.frame(ef)

ef$IDEOL_f <- factor(ef$IDEOL_f, 
                       levels = c("right", "center", "left"), 
                       labels = c("right", "center", "left"))

colors = c("black", "grey50")
m2_pred <- ggplot(ef, aes(IDEOL_f, fit, color=name_f, fill=name_f)) + geom_point(size=2) + 
  geom_errorbar(aes(ymin=lower, ymax=upper), width=0.3, size=1) + 
  theme_bw(base_size=12) +
  theme_bw() + labs(x="", y="Probability of preferring candidate", color="Candidate name", title="B") +
  scale_colour_manual(values = colors, aesthetics = c("colour", "fill")) + 
  guides(fill=FALSE) + coord_flip() +
  scale_x_discrete(limits = rev(levels(ef$IDEOL_f))) +
  theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
  theme(plot.title = element_text(size=11)); m2_pred

### Marginal effects

m2_marg <- summary(margins(m2, at = list(IDEOL_f = c("left","center", "right")), variables = "name_f", data=d_long))

colors_2 <- c("black")

m2_marg$Ideology <- factor(m2_marg$IDEOL_f , 
                            levels = c("left", "center", "right"), 
                            labels = c("left", "center", "right"))

m2_eff <- ggplot(m2_marg, aes(x = Ideology, ymin=lower, ymax = upper,
                              color=factor)) + geom_errorbar(color=c("black"), width=0.3, size=1, 
                                                             position=position_dodge(.5)) + 
  geom_point(aes(x=Ideology, y=AME, color=factor), color=colors_2, size=2) +
  theme(legend.position="none") +
  coord_flip() +
  labs(x="", y="Effect of non-Swiss name on probability of being preferred", title="A", color = "Candidate name") + 
  theme_bw() +
  scale_x_discrete(limits = rev(levels(m2_marg$exp_arm_f))) + 
  geom_hline(yintercept=0, linetype="dashed", size=1)  +
  scale_fill_manual(values=colors_2) + theme(legend.position="right"); m2_eff

m2_eff_pred <- grid.arrange(m2_eff, m2_pred, 
                            ncol=2, nrow = 1,
                            widths = c(1.8, 2.2)); m2_eff_pred 


# ---------------- All candidate characteristics

# ----- Z0

margins_all <- margins::margins(m3, at = list(exp_arm_f = c("0", "1")))
margins_all_d <- summary(margins_all)

margins_all_0_d <- margins_all_d %>% filter(exp_arm_f=="0")

modelFrame <- data.frame(Variable = margins_all_0_d$factor,
                         Coefficient = margins_all_0_d$AME,
                         SE = margins_all_0_d$SE)

modelFrame <- modelFrame[-4, ]

# Add reference categories

Variable <- c("age_35", "policy_jugendpolitik", "profession_hausarzt", 
              "children_2", "name_moser")
Coefficient <- c("", "", "", "", "")
SE <- c("", "", "", "", "")

ref_cat.data <- data.frame(Variable, Coefficient, SE)

allModelFrame <- merge(modelFrame, ref_cat.data, all=TRUE)

#Numeric variables

allModelFrame <- allModelFrame %>% map_df(rev)

allModelFrame$Coefficient <- as.numeric(as.character(allModelFrame$Coefficient))
allModelFrame$SE <- as.numeric(as.character(allModelFrame$SE))

allModelFrame$Variable <- factor(allModelFrame$Variable, levels=order_var)

#Coefficient plot
coef_conjoint_z0 <- PlotCoefficients(z)
coef_conjoint_z0 <- PlotCoefficients(z) + scale_y_continuous("Probability of preferring candidate", limits = c(-0.3, 0.3)) + scale_x_discrete("", labels = my_labels) + labs(title="No competence information") +
  my_theme + theme(axis.text.y = element_text(face = c("plain", "bold", "plain", "plain", "plain", "bold", "plain", "plain", "bold", "plain", "plain", "plain", "plain", "plain", "bold", "plain", "bold"))) +
  scale_colour_manual("Model", values = "black") + theme(legend.title = element_blank())  +
  theme(plot.title = element_text(size=11))
coef_conjoint_z0

# ----- Z1

margins_all <- margins::margins(m3, at = list(exp_arm_f = c("0", "1")))
margins_all_d <- summary(margins_all)

margins_all_1_d <- margins_all_d %>% filter(exp_arm_f=="1")

modelFrame <- data.frame(Variable = margins_all_1_d$factor,
                         Coefficient = margins_all_1_d$AME,
                         SE = margins_all_1_d$SE)

modelFrame <- modelFrame[-4, ]

allModelFrame <- merge(modelFrame, ref_cat.data, all=TRUE)

#Numeric variables

allModelFrame <- allModelFrame %>% map_df(rev)

allModelFrame$Coefficient <- as.numeric(allModelFrame$Coefficient)
allModelFrame$SE <- as.numeric(allModelFrame$SE)

allModelFrame$Variable <- factor(allModelFrame$Variable, levels=order_var)

#Coefficient plot
coef_conjoint_z1 <- PlotCoefficients(z)
coef_conjoint_z1 <- PlotCoefficients(z) + scale_y_continuous("Probability of preferring candidate", limits = c(-0.3, 0.3)) + scale_x_discrete("", labels = my_labels) + labs(title="Competence information") +
  my_theme + theme(axis.text.y = element_text(face = c("plain", "bold", "plain", "plain", "plain", "bold", "plain", "plain", "bold", "plain", "plain", "plain", "plain", "plain", "bold", "plain", "bold"))) +
  scale_colour_manual("Model", values = "black") + theme(legend.title = element_blank()) +
  theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
  theme(plot.title = element_text(size=11))
coef_conjoint_z1

# ----- Difference Z1-Z0

modelFrame <- data.frame(Variable = rownames(summary(m3)$coef),
                         Coefficient = summary(m3)$coef[, 1],
                         SE = summary(m3)$coef[, 2])

modelFrame <- modelFrame[-c(1:14), ]

allModelFrame <- merge(modelFrame, ref_cat.data, all=TRUE)

#Numeric variables

allModelFrame <- allModelFrame %>% map_df(rev)

allModelFrame$Coefficient <- as.numeric(allModelFrame$Coefficient)
allModelFrame$SE <- as.numeric(allModelFrame$SE)

order_var <- c("exp_arm_f1:children_d3","children_2", 
               "exp_arm_f1:policy_fWirtschaftspolitik", "exp_arm_f1:policy_fUmweltpolitik", "exp_arm_f1:policy_fSicherheitspolitik", "policy_jugendpolitik", 
               "exp_arm_f1:age_f61", "exp_arm_f1:age_f48", "age_35", 
               "exp_arm_f1:profession_fVerkäufer", "exp_arm_f1:profession_fKaufmännischer Angestellter",
               "exp_arm_f1:profession_fSozialarbeiter", "exp_arm_f1:profession_fPrimarlehrer", 
               "exp_arm_f1:profession_fRechtsanwalt", "profession_hausarzt", 
               "exp_arm_f1:name_fNon-Swiss name", "name_moser")

allModelFrame$Variable <- factor(allModelFrame$Variable, levels=order_var)

#Coefficient plot
coef_conjoint_z01 <- PlotCoefficients(z)
coef_conjoint_z01 <- PlotCoefficients(z) + scale_y_continuous("Probability of preferring candidate", limits = c(-0.3, 0.3)) + scale_x_discrete("") + labs(title="Difference") +
  my_theme + theme(axis.text.y = element_text(face = c("plain", "bold", "plain", "plain", "plain", "bold", "plain", "plain", "bold", "plain", "plain", "plain", "plain", "plain", "bold", "plain", "bold"))) +
  scale_colour_manual("Model", values = "black") + theme(legend.title = element_blank()) +
  theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
  theme(plot.title = element_text(size=11))
coef_conjoint_z01

library(gridExtra)
coef_conjoint_z <- grid.arrange(coef_conjoint_z0, coef_conjoint_z1, coef_conjoint_z01,
                            ncol=3, nrow = 1,
                            widths = c(2.2, 1, 1))
coef_conjoint_z



